gusucode.com > 支持向量机工具箱 - LIBSVM OSU_SVM LS_SVM源码程序 > 支持向量机工具箱 - LIBSVM OSU_SVM LS_SVM\stprtool\svm\kernelsk.m
function [Alpha,bias,sol,t,flps,margin,up,lo,trn_err]=...
kernelsk(X,I,epsilon,ker,arg,tmax,C,up,lo)
% KERNELSK kernel Schlesinger-Kozinec's algorithm, SVM (L2).
% [Alpha,bias]=kernelsk(X,I,epsilon,ker,arg,tmax,C)
%
% [Alpha,bias,sol,t,flps,margin,up,lo,trm_err]=...
% kernelsk(X,I,epsilon,ker,arg,tmax,C)
%
% This algorithm solves the SVM problem with the quadratic
% penalizition of classification violations using the
% kernel Schlesinger-Kozinec's algorithm.
%
% Inputs:
% X [NxL] training patterns, N is dimension and L number of patterns.
% I [1xL] labels, 1 for 1st class and 2 for 2nd class.
% epsilon [1x1] precision of found solution. The margin of found
% hyperplane is less than the optimal margin at most by epsilon.
% ker [string] kernel, see 'help kernel'.
% arg [...] argument of given kernel, see 'help kernel'.
% tmax [int] maximal number of iterations.
% C [real] trade-off between margin and training error.
%
% Outputs:
% Alpha [1xL] Weights (Lagrangians) of patterns.
% bias [real] bias (threshold) of found decision rule.
% sol [int] 1 solution is found
% 0 algorithm stoped (t == tmax) before converged.
% -1 hyperplane with margin greater then epsilon
% does not exist.
% t [int] number of iterations.
% margin [real] margin between classes.
% flps [int] number of used floating point operations.
% up [1,t] evolution of the upper bound on the optimal margin.
% lo [1,t] evolution of the lower bound on the optimal margin.
% trn_err [real] training error (empirical risk).
%
% See also KERNELSKF, SVM.
%
% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 02.11.1999, 13.4.2000
% Modifications
% 23-October-2001, V.Franc
% 19-Septemberr-2001, V.Franc, comments changed.
% 18-August-2001, V.Franc, m^* - m <= epsilon instead of <= epsilon/2
% added uppper and lower bound on
% 17-August-2001, V.Franc, renamed to KERNELSK.
% 13-July-2001, V.Franc, comments
% 12-July-2001, V.Franc, C, bias and normal vect. normalized.
% 11-July-2001, V.Franc, Rosta Horcik proved that the computation
% of threshold is OK.
% 10-July-2001, V.Franc, derived from kekozinec2
flops(0);
% set default values of the input argiments
if nargin < 7,
C = inf;
end
% maximal number of iteraions
if nargin < 6,
tmax = inf;
end
% indexes of pattens in the 1st and 2nd class
xinx1 = find(I == 1);
xinx2 = find(I == 2);
X1=X(:,xinx1); % patters from 1st class
X2=X(:,xinx2); % patters from 1st class
l1 = size(X1,2); % number os patterns
l2 = size(X2,2);
% compute kernel matrices
K1 = kernel( X1, X1, ker, arg ); % [l1 x l1]
K2 = kernel( X2, X2, ker, arg ); % [l2 x l2]
K12 = kernel( X1, X2, ker, arg ); % [l1 x l2]
% make problem lin-separable in high dimensional space
if C ~= 0,
CD1=eye(l1,l1)/(2*C); % additional diagonal matrix
CD2=eye(l2,l2)/(2*C);
K1=K1+CD1;
K2=K2+CD2;
end
% convex coeficients defining normal of the decision hyperplane
% (they correspond to the Lagrangian multiplyers).
s1 = zeros(l1, 1);
s2 = zeros(l2, 1);
% initial solution
s1(1)=1; % take the 1st pattern from the 1st class
s2(1)=1; % take the 2nd pattern from the 2st class
sol=0;
t = 0;
up=[];
lo=[];
% main cycle
while sol == 0 & tmax > t,
t = t + 1;
sol = 1;
% -- compute auxciliary variables --
a = s1'*K1*s1;
b = s2'*K2*s2;
c = s1'*K12*s2;
f = K2*s2;
e = K1*s1;
d1 = e - K12*s2;
d2 = f - (s1'*K12)';
[d1min,inx1] = min(d1);
[d2min,inx2] = min(d2);
% projection x \in X_1 on (w_1 - w_2)
proj1 = (d1min + b -c )/sqrt(a-2*c+b);
% projection x \in X_2 on (w_2 - w_1)
proj2 = (d2min + a - c)/sqrt(a-2*c+b);
if sqrt( a -2*c +b) <= 0,
% algorithm has converged to the zero vector --> classes overlap
sol = -1;
break;
end
% --- compute stop condition for the alpha1 (1st class) -----
% (proj1 < proj2) ~ the worst point will be used for update
if (proj1 < proj2) & (proj1 <= (sqrt(a-2*c+b) - epsilon)),
% -- Adaptation phase of vector alpha1 ----------------------------
k = (a - d1min - c)/(a + K1(inx1,inx1) - 2*e(inx1) );
k = min( 1, k );
s1 = s1 * (1-k);
s1(inx1) = s1(inx1) + k;
sol = 0;
else
% --- compute stop condition for the alpha2 (2st class) ------
if proj2 <= (sqrt(a-2*c+b) - epsilon ),
% -- Adaptation phase ----------------------------------
k = (b - d2min -c)/(b + K2(inx2,inx2) - 2*f(inx2) );
k = min( 1, k );
s2 = s2 * (1-k);
s2(inx2) = s2(inx2) + k;
sol = 0;
end
end
% store up=||w||/2 and current margin m(w1-w2,theta) = min( m1, m2)
m = min([proj1,proj2]) - 0.5*sqrt(a-2*c+b);
up = [up,sqrt(a-2*c+b)/2];
lo = [lo,m ];
end
% does the found hyperplane separate the classes or not ?
if sol == 1 & ((d2min + a - c) < 0 | (d1min + b -c ) < 0),
sol = -2;
end
% --- determine threshold --------
%theta = 0.5*( s1'*K1*s1 - s2'*K2*s2 );
% sqared margin in transfromed space
margin2 = s1'*K1*s1 - 2*s1'*K12*s2 + s2'*K2*s2;
% threshold after normalization
theta = ( s1'*K1*s1 - s2'*K2*s2 )/margin2;
% solution (normal vect. in the transformed space) after normalization
s1=2*s1/margin2;
s2=2*s2/margin2;
% --- compute margin and classify training patterns
if C~=0,
K1=K1-CD1;
K2=K2-CD2;
end
margin = 1/sqrt( s1'*K1*s1 - 2*s1'*K12*s2 + s2'*K2*s2 );
dpred1 = (K1*s1 - K12*s2 )' - theta;
dpred2 = (K12'*s1 -K2*s2 )' - theta;
% classification error on the traning set
trn_err = length( find([dpred1,-dpred2] < 0) )/(l1+l2);
%----------------------
% make SVM-like output
Alpha=zeros(1,l1+l2);
Alpha(xinx1)=s1;
Alpha(xinx2)=s2;
bias = -theta;
% overall number of used float point operations
flps=flops;
return;